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We present a flexible method that can calculate Bloch modes, complex band structures, and impedances of 
two-dimensional photonic crystals from scattering data produced by widely available numerical tools. The 
method generalizes previous work which relied on specialized multipole and FEM techniques underpinning 
transfer matrix methods. We describe the numerical technique for mode extraction, and apply it to calculate a 
complex band structure and to design two photonic crystal antireflection coatings. We do this for frequencies 
at which other methods fail, but which nevertheless are of significant practical interest. 



I. INTRODUCTION 

When modeling photonic crystals (PCs), it is impor- 
tant to consider all the relevant Bloch modes. Light at 
a fixed frequency, polarization, and incident angle exists 
in a PC as a superposition of a set of propagating and 
evanescent Bloch modes, the PC's eigenstates. At low 
frequencies, only one mode generally needs to be con- 
sidered. For light at frequencies above the first Wood 
anoma Ijiil, each row of holes in the PC diffracts light into 
several propagating orders, so the PC may support mul- 
tiple propagating Bloch modes. At the PC's front and 
back interfaces, some of its modes couple via reflection, 
affecting the overall reflection and transmission through 
the PC, so it is important to model all relevant modes. 

It is often important to include evanescent modes^. If 
the PC is not long — for example, if it is a layer in a thin 
antireflection coating — then evanescent modes can play 
a role in energy transporlP. Evanescent modes can also 
play a role in field matching across an interface between 
PCP or PC waveg uideP. The propagative qualities of 
an evanescent mode are well-represented by its complex 
band structur^, which augments the traditional band 
structure, conveying information about the rate at which 
the mode accumulates phase together with information 
about the mode's decay rate. 

There have been a number of studies seeking to de- 
rive impedance-like quantities to characterize reflection 
at PC interfaces by a scalail^SI. Furthermore, a number of 
studies have adapted metamaterial parameter extraction 
technique^ to photon ic cry stals, and used them to design 
antireflection coating^iSIll], However, since these tech- 
niques characterize reflection and transmission by a sin- 
gle complex number each, they cannot handle problems 
involving multiple modes, where every mode reflects into 
every other mode. Scalar-based methods generally give 
manifestly incorrect results for light at frequencies above 



the first Wood anomaly, which ranges from a^/X = 1/n 
for normally incident light to a^/A = l/2n for light at the 
Brillouin-zone edge, where is the length of the lattice 
vector parallel to the interface, A is the free space wave- 
length and n is the PC's background index. Above this 
frequency, generally several Bloch modes must be simul- 
taneously considered in each PC, regardless of whether 
these modes are propagating or evanescent. Reflection 
at a PC/PC interface is well-described by a matrix that 
maps incident mode s to refiected modes, as we have 
shown previousljffl^. In our experience, the minimum 
acceptable dimension of this reflection matrix, as argued 
in Sec. |II A[ is usually 
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where 9i is the incident angle from a uniform dielectric 
with the PC's background index, and [a:J denotes the 
floor of X. We have previously achieved accurate results 
modeling PC stack s using impedance matrices of this di- 
mension and higheif^li^ESI. 

A number of methods for finding multiple Bloch modes 
and complex band structures have been demonstrated. 
Transfer-matrijP^ and scattering-matri>l^ based meth- 
ods were developed to derive a PC's Bloch modes from 
the properties of a single grating layer. The plane wave 
expansion method has also been extended to include 
evanescent modeJ^^. Finally, Ha et al. presented a 
method for extracting Bloch modes from the o utput of an 
EM solvei^i^, or even near-field measurement^iSMI 
improve the accuracy, stability and efficiency of Ha et 
aVs method and extend it to calculate PC impedances 
for two-dimensional (2D) PCs, which can be used to cal- 
culate reflection and transmission at interface d^ l ^^ l These 
PC impedances and the reflection and transmission oper- 
ators are represented by matrices; our method supports 
the presence and interaction of multiple Bloch modes and 
so it can work well both above and below the flrst Wood 
anomaly. 

We have made software available that uses the method 
described in this paper to calculate PCs' Bloch modes. 
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complex band structures, and impedances. The soft- 
ware, called BlochCode, can then use these complex band 
structures and impedances to calculate reflection and 
transmission matrices and coefficients for arbitrary stacks 
of PCs. BlochCode is open-source and is available on the 
internetP. 

In Sec. |lT] we present our method for finding Bloch 
modes from the electric field E and the magnetic field H 
in a PC structure. Sec. |II A| recaps some useful results 



from our previous worlf^^ and provides some background 
theory. Sec. |II B| details our improvements to Ha et aZ.'s 
mcthod^^ of finding Bloch factors and modal fields, and 



Sec. Ill outlines our procedure for successfully applying 
this method to minimize the residual derived in Sec. Ill Bl 
Sec. |IIC| explains how we calculate PC impedance matri- 
ces from the modal fields. In Sec. IV we apply our method 



to demonstrate its utility. In Sec. IV A| we calculate the 
complex band structure for light normally incident on a 
triangular lattice PC. In Sec. |IVB] we reproduce the de- 
sign process of a known antireflection coating for a PC, 
at a frequency and incident angle for which it is critical 
to include at least two Bloch modes in the calculations. 
Finally, in Sec. |IV C| we use our method to design an 
all-polarization antireflection coating for a square lattice 
self-collimating PC, at a high frequency where a scalar 
method cannot find a coating for the PCl^il. 



II. THEORY 

Our method uses a two-step process to extract a PC's 
modes and impedance from the field in a finite length of 
the PC. The PC is assumed to be two-dimensional, loss- 
less, and to have relative permeability — 1. Like Ha 
et al.'s methocP^, we could use data generated by FEM 
or FDTD simulations, or even experimentally measured 
by a near-field probe such as a SNOMiSi^ although the 
impedance part of our method is not valid for SNOM 
data, which is derived from a 3D object. First, the Bloch 
factors and the Bloch modal fields are found (Sec. II B I, 



then these modes are analyzed to calculate the PC's 



impedance (Sec. II C) 



A. Background Theory 

Two-dimensional PCs in the x — y plane may be de- 
scribed as a stack of gratings parallel to the x axiJ^H^ 
each of which diffracts incident light into an infinite set 
of grating orders. At the edge of each unit cell, the PC's 
Bloch modes may be written as a superposition of the 
underlying grating orderJ^^. Their directions are given 
by the grating equation 
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where is the x component of the incident plane wave's 

(p) 

wavevector, fci is that of the pth diffraction order, and 



is the length of the lattice vector parallel to the x-axis. 
The wavevector component in the direction perpendicu- 



lar to the grating is k. 



(p) 
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where k is the 



wavenumber in the medium. Evanescent grating orders 

(p) • (p) 

have imaginary ky , so for a given k and kx , the number 

of propagating grating orders is the number of solutions 
to Eq. ^ with real or Af^in in Eq. In our 

experience, Mmin also provides an upper bound on the 
number of propagating Bloch modes, and at non-normal 
incidence is a lower bound on the number of Bloch modes 
required to model a PC accurately. At normal incidence, 
symmetry allows odd modes to be ignored, so in this 
case good results may be obtained with fewer than M^un 
modes — see Sec. |IVC| Using Bloch modes found from ac- 
curate multipole and FEM transfer matrix mcthods^"^ 
we have consistently had success modeling PCs with no 
more than M,„in -t- 2 Bloch modes. 

Bloch's theorem relates the electric and magnetic fields 
associated with each mode at equivalent points in differ- 
ent unit cells of a PC. The ratio of each mode's field 
at points separated by the lattice vector ei — {ax,0) 
is e*'^^'^^. For the PC's other lattice vector 62, this ra- 
tio is different for each mode and is the mode's Bloch 
factor, denoted by /i. Calculating /i for each mode is 
the goal of Sec. |II B[ For square and rectangular lattices, 
62 = (0, Cj,) and /i — e*'"""", where ky is the y compo- 
nent of the mode's wavevector. For triangular lattices, 
the lattice vector 62 is (aa;/2, ay) and so the Bloch factor 
may be written /i = e 



i{k^a^/2+kyay) 



Bloch modes come in forward/backward pairs. Popov 
et al. provide a useful discussion of symmetry 
properties'-^. We assume mirror symmetry in each unit 
cell, which means that each backward mode's field profile 
in a unit cell is the reflection on the x-axis of its forward 
partner's. The Bloch factors of a pair are related because 
of this: for square and rectangular lattices, /if, = 
where fif and /Xb are respectively the Bloch factors of 
the forward and backward modes. For triangular-like 
lattices, the symmetry is more complicated since the re- 
flection of 62 is not —62, the translation corresponding 
to the field ratio but {ax/2,—ay); these vectors 

differ by — ei. Accounting for this discrepancy, we find 
/ib = e"**^^"^//!/ for triangular lattices. 

A PC's impedance is defined in terms of two matrices, 
E and Yot E = polarized light, each matrix 

maps a vector of forward Bloch mode amplitudes c+ to a 
vector of the Ez or fields associated with each grating 
diffraction order. Specifically, -E^.m, the (p, TO)th element 
of E, is the E^ field of normalized mode m due to forward 
and backward plane waves in grating order p, at the cen- 
tre {x = 0) of a unit cell's edge. Thus, for a set of forward 
propagating/decaying Bloch modes c+, the field compo- 
nents along the edge of the unit cell, i.e., the quantities 
that are continuous across an interface between PCs or 
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dielectrics, are 



(3) 

where Ep and Hp are the rows of E and H corresponding 
to grating order p. In the H = polarization, E and 
H map to Ex and Hz fields, and these quantities replace 
Ez and in Eq. (|3|. 

Previousljff^, we defined PC impedances in terms of 
these matrices. For Ez polarized light, the impedance of 
a PC is 

Z = Ho^(I + Q)E + Eo^(I-Q)H, (4) 
and for Hz polarized light it is 

Z = -(Ho^(I-Q)E + Eo^(I + Q)h), (5) 



where E and H are calculated for the PC, and Eq and 
Ho are calculated for a reference material, usually free 
space. Q is a diagonal matrix that takes into account 
the half-period shift of gratings in triangular lattice PCs: 
for square lattices Q = I, and for triangular lattices Q = 
diag((— where p is the grating order. 

Given impedances Zi and Z2 for two PCs, it is sim- 
ple to calculate the reflection and transmission matrices 
across their interfac^^^ 
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where A12 



B. Finding modes 

Our method of finding the Bloch modes and Bloch fac- 
tors is based on the method presented by Ha et al}'' , al- 
though our method offers some significant improvements 
in accuracy and efficiency. We take field data for several 
unit cells of a PC, and try to write it as a superposition of 
Bloch modes, thus finding the modal fields and Bloch fac- 
tors. The final steps of our mode-finding method impose 
symmetry relationships between forward and backward 
modal fields, increasing accuracy by almost halving the 
number of unknowns in the problem. We now outline our 
method. 

In an EM solver, we simulate a section of 2D PC 
with Bloch-Floquet periodic boundary conditions on two 
boundaries, and uniform dielectric on the others (Fig.jl]). 
We sample the Ez or E^ (depending on polarization) field 
component at many {Np) points in unit cell ^ = 0, and 
then at the equivalent points in each of the other unit 
cells. If desired, Ey, H^, Hy, or Hz may be used in 
place of or in addition to Ez and E^- For triangular 
lattice PCs, we use the field in the simulated unit cells 
(dashed edges in Fig. [T]) to calculate the field in the unit 
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FIG. I: Schematic of L = 5 PC structures for a square 
and a triangular PC lattice. The squares with solid 
edges are the unit cells used by our method. For the 
triangular lattice PC, the field in the solid-edge unit 

cells are calculated from the unit cells of the simulated 
structure (dashed edges) using Bloch's theorem, with 
the ratio 6'*^="°^ between adjacent cells' fields. 



cells separated by a lattice vector (solid edges); we ap- 
ply Bloch's theorem with integer multiples of the lattice 
vector (Oa:, 0). 

We seek to write these electric field components as a 
superposition of forward and backward Bloch modes. So 
we want to express every Ui{y), i.e., the Ez or E^ field 
component for sampled point r in unit cell £, as 



(l^„-^-i-^)A„'(r)-Hu;(£,r), 

(7) 

where ^m(r) and /i™ are respectively the modal field 



and the Bloch factor of forward mode 



denotes 



backward modes, and w(£, r) is the residual error. More 
specifically, for forward modes, Am{v) is the field compo- 
nent of mode m at point r of the first unit cell, € = 0. The 
Bloch factor /i^ is the ratio of the field in cells l+l and 
so /if„Am(r) is the field component of forward mode 
m at point r of unit cell £. To avoid ill-conditioning, the 
field Am' (r) at point r of each backward mode m! is de- 
fined in the last unit cell, i = L—1. This means that the 
coefficients of ^^(r) and ^m'(r) in Eq. (TJ) have moduli 
no greater than 1. As noted in Sec. |II A the Bloch fac- 
tor ^rri' of each backward mode is related to that of its 
forward partner; we enforce this relationship in practice, 
thereby halving the number of Bloch factors that must 
be found. 

Equation ([t]) for all I and all sampled r may be written 
in matrix form as: 



U ^ CA + W, 



(8) 



where U contains the Ez or E^ field components from the 
EM solver, A is a matrix of modal fields, C is a matrix 
constructed from Bloch factors, and W is a matrix of 
residuals w{l,r) that must be minimized. U is a L x 
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Np matrix: the field in its £th row and rth column is 
Ue^r = Ui{r), the field component at point r in unit cell 
£. Similarly, A is a M x Np matrix; the field in its mth 
row and rth column is Am,r — ^)n(r), the field of mode m 
at point r in cell ^ = for forward modes, or cell £ — L — 1 
for backward modes. C is a L x Af matrix. For a forward 
mode m, the {£,m)th element of C is Hm^, and for a 

backward mode m', the (£, rn')th element is ^ ^• 

If multiple field components (e.g. E^, and Hy) are to 
be used to find the modes, then the additional data can 
be added as extra columns in U. 

We start the optimization process knowing U, and 
with information about the structure of C, and no di- 
rect information about A. In our method, we first find 
the Bloch factors that determine C, a relatively difficult 
problem. Once C is known, solving Eq. (|8| for the modal 
fields A becomes a pure least-squares problem that can 
be solved accurately and efficiently using standard tech- 
niques. 

To find the modes, we seek to minimize the difference 
between the observed field U and the superposition of 
Bloch mode fields CA. That is, we seek to minimize 
||W||| in Eq. the sum of squared moduli of the el- 
ements of W. Constraining the problem by dividing by 
the squared Frobenius norm | |U| \p of U, the quantity we 
minimize is 



improve the quality of our results by enforcing this rela- 
tionship in the minimization process. 

We commence by partitioning the forward (/) and 
backward (6) modes, and the points in the left (L; 
y < ay/2) and right {R; y > ay/2) halves of the unit 
cell: 



U = (Ui,Ufl), C = (C;,C 



A = 



Al,6 



(12a) 



(12b) 



After normalization, the field of a backward mode is the 
field of its forward partner reflected about the x-axis, 
thus 



(AL,h, A 



R.b) 



7AfljP,7A 



(13) 



where P is the permutation matrix that maps points 
y) to {x, y), and 7 is a normahzing diagonal ma- 
trix whose elements are the ratio of backward and for- 
ward mode amplitudes. The columns of Arj and A^j^f,, 
corresponding to points in the right half of the unit cell, 
can easily be ordered so that P = I; from now on we 
assume this ordering. Eq. Q can now be written with 
roughly half as many unknowns. 
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U|||,. First we eliminate A from 



where w"^ 

Eq. ([9| in order to find C with a numerical minimizer 
We use an alternative representation of the Frobenius 
norm, ||U||i? ~ ■\/tr(TP^Tj), to write 



tr((U^ - A^C^)(U - CA)) 



luill^ 



(10) 



Finding A for arbitrary C is a standard least-squares 
problem; the optimal A satisfies C^CA = C^U. We 
expand Eq. ( 10 ), twice apply this relation, and rearrange 
to get 



2 _ tr(U^CC+U) 



(11) 



where C"*" = (C^C)~^C^ is the Moore-Penrose pseu- 
doinverse of C. 

Using Eq. (11) and a numerical minimizer, the Bloch 
factors that determine C may often be found to a useful 



level of accuracy (see Sec. Ill for implementation details) 



In order to improve the accuracy and reliability of the 
results, we impose further ph ysica l constraints. 

The PC impedance methocPE^l assumes the unit cell 
to be up-down symmetric, which causes the forward and 
backward modes to be related. So far, we have only im- 
posed a relationship between the forward and backward 
Bloch factors, not the modal fields within each unit cell. 
We can halve the number of unknowns in A and strongly 



(Ui,UH) = (C/,C,7)U^;-^ A 



Al,/ Arj 



w. 



(14) 



Cb7 represents each backward mode's amplitude in each 
cell, relative to that of the corresponding forward mode 
in cell 0. 



The constraints on A (Eq. (13)) mean that Eq. (14) 



does not have a least-squares form, so may not be imme- 
diately simplified in the way that Eq. ^ led to Eq. ( [Tl] ) . 
To transform Eq. ( 14 ) into a more useful form, we block- 
diagonalize A and right-multiply by the matrix ( i ) , 
to show 



(U+, U_) = (C+A+, C_A_ 



w. 



(15) 



Here we have introduced the symmetric and antisym- 
metric forms U± = Ul ± Vr, C± = C/ ± 0^7, and 
A± = Alj ± Arj. 



Eq. ( 15 ) takes the form of two independent least- 



squares equations, each with half the dimension of 
Eq. (14 1. The two equations must be satisfied simul- 



taneously, so to find the Bloch factors we can minimize 
2 ||U+-C+A+||| + ||U_-C_A 



|2 

If 



|U+ 



lU IP 



or equivalently 

tr(Uf C+CjU+) + tr(U^C_Cj:U_) 



1 



|u+|||. 



(16) 



(17) 



Again, this quantity may be minimized by a numerical 
optimizer. The residual w'^ for any solution to Eq. (17) 
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is equal to the residual obtained by inserting the solution 
into Eq. (Ill: the two equations differ only in the symme- 



try constraint on backward modal fields (Eq. ( 13 )). Com- 
pared to Eq. ( 11 1, we have removed NpM unknowns from 
A (where Np ^ M is the number of sampled points in 
each unit cell) , halving its dimension at the cost of adding 
M unknowns to C± as 7. These new unknowns must be 
found simultaneously with the Bloch factors using a nu- 
merical minimizer, so it is important to supply a good 
starting estimate; our method for doing so is detailed in 

secnni 



C. Calculating impedance 

Once the Bloch factors and 7 are known, the modal 
fields can be reconstructed and analyzed to determine 
the PC's impedance. The essential quantities for this cal- 
culation are the E and H field components in the plane 
of the PC interface (i.e., Ez and H^, or E^ and H^, de- 
pending on polarization) of each Bloch mode m along the 
left edge {y = 0) of a unit cell (see Fig. [l]). These quan- 
tities, Em(x) and Hm{x), may be found from Eq. (15) 
using the known values for C+ and C_ and inserting the 
appropriate E 01 H fields into U+ and U_. 



To calculate the impedance, we find the E and H ma- 
trices for the PC, as defined in Se c. |II A Inserting mul- 
tiples of unit vectors c_|_ into Eq. (|3|), we can show that 



Emi^x) — •^m, ^ ^ Ep 



ifc<P'a 



(18a) 



(18b) 



where Am is the amplitude of the normalized mode m, 
and Ep^m and Hp „i are the elements of E and H. It is 
straightforward to exploit the orthogonality of the plane 
wave grating diffraction orders to show that 



ax/2 



AmEp^m^l/a^ I Err,{x)e '^-^^ dx, (19a) 

-Ox/2 



ax/2 



AmHp.m^l/a^l ff™(x)e-*'=^''^da;. (19b) 

J-ax/2 



Eqs. ( 19 ) let us calculate each element of the E and H 
matrices, up to a normalization constant A„l per col- 
umn. We remove the constants by calculating the PC's 
impedance (Eq. (|4| or ([5])) with the PC itself as the refer- 
ence material: by reciprocity-derived Bloch mode orthog- 
onality relational^, this quantity should be the identity 
matrix. The diagonal entries of this matrix are the Am^] 
the off-diagonal terms, which should be zero, provide an 
error estimate. After normalizing the E and H matrices 
for the PC, we calculate its impedance matrix Z from Eq. 
(|4]) or ([5]) using a reference medium such as free space. 



III. NUMERICAL PROCEDURE 

Having outlined the theoretical basis of our method 
for finding the Bloch factors and impedance of a PC at a 
given frequency, incident angle, and polarization, we now 
provide some practical detail about our implementation 
of the method. We outline the procedure for Af = 3 pairs 
of Bloch modes. 

In COMSOL Multiphysics 4.2, we simulate a 1 x 8 unit 
cell sample of PC, embedded in its background dielectric, 
with Bloch-Floquet periodic boundary conditions along 
the two long boundaries (Fig. [l] shows a 1 x 5 struc- 
ture). Eq. (15) is a set of LNp equations, with 2M 
and MNp unknowns in C± and A± respectively. To be 
overspecified, the method requires LNp > MNp + 2AI; 
thus L = 8 periods and a large Np is sufficient to find 
M — 3 modes. A deeper structure with more unit cells 
does not necessarily provide useful information about ad- 
ditional evanescent modes, as their amplitude deep inside 
the structure may be negligible. From COMSOL we ex- 
port the relevant E and H field components in the L = 8 
unit cells, sampled over a 101 x (50L + 1) grid. 

In order to compute a mode, it must be present in the 
structure with sufficient amplitude to be detected. Light 
at normal incidence often fails to excite odd Bloch modes; 
these uncoupled mode^^ consequently cannot be found 
by an optimization, which loses accuracy in searching for 
modes that are not present. At frequencies above the first 
Wood anomaly, the frequencies at which the higher order 
modes are most important, this problem may be avoided 
by exciting the PC slab not with a normally incident 
plane wave, but with the first grating diffraction order. 
This technique is used in Sec. |IV A| and Sec. |IV C[ If the 
uncoupled mode is not relevant to a particular problem, 
it may instead be ignored. 

If we seek to find M — 'i Bloch modes, then find- 
ing a global minimum of Eq. (17) involves searching 
for 2M = 6 complex numbers. This is a hard prob- 
lem if attacked directly, but we use an algorithm that 
gives more consistent success by providing a good start- 
ing estimate. We start by minimizing the residual 



in Eq. (11), which forces a relationship between forward 



and backward Bloch factors but not the modal fields. 
This involves finding only M complex numbers. As a 
starting estimate for the forward Bloch factors, we either 
take the result of a neighboring simulation, or the ana- 
lytically calculated Bloch factors for the dielectric back- 
ground of the PC. At every step of the minimization, 
evanescent modes are sorted into forward and backward 
decaying modes, based on the moduli of their Bloch fac- 
tors. The minimization can be done by any standard 
numerical minimizer, such as SciPy'^ fmin, which is a 
modified Nelder-Mead optimizatiorPSl. At this point, the 
results are equivalent to those from the method of Ha et 
aZ.P^^, except that we have lessened the likelihood of C be- 
ing ill-conditioned by renormalizing the backward Bloch 
factors jjirn' in Eq. ([7| and setting their phase origin to 
the end of the PC. 



6 



Occasionally, we encounter an instability in which a 
pair of modes have very large equal and opposite field 
amplitudes and very small Bloch factors. When this oc- 
curs, we follow a Gram-Schmidt-like process: we sub- 
tract the field of non-problematic modes (i.e., modes with 
> 10^'^) from U and repeatedly minimize Eq. (11 1 to 
find each of the remaining modes individually. 



Using the solution to Eq. (11) as our estimate for the 



Bloch factors, the modal fields may be found with a least- 
squares optimization. The average field ratio of each pair 
of backward and forward modes gives us an estimate for 
7. We now have a plausible estimate for 7 and the Bloch 
factors, which we can use as a starting estimate to mini- 
mize Eq. (fTTl). 



To further refine the estimates, we repeatedly iterate 
through the modes, fixing all but one fi and the corre- 
sponding element of 7, minimizing Eq. (17 1 to find the 
two variables. After this process, we finally minimize 



Eq. (17) across all 6 complex dimensions simultaneously 
to obtain the correct Bloch factors and modal fields from 
which we calculate impedances. Forward and backward 
propagating modes are sorted based on their flux^^, be- 
fore impedances are calculated as outlined in Sec. |II C| 



IV. APPLICATIONS 

We now apply our method to a range of typical prob- 
lems. Each of these problems involves frequencies above 
the first or second Wood anomaly — frequencies at which 
scalar methods fail and multiple modes are required to 
describe the system. BlochCode, software that imple- 
ments our method in Python, using SciPj^'^ and bage , 
is freely available on the internetPSl we use it here. 



A. Complex band structure 

The first application of our method is to calculate the 
complex band structure of a PC. The PC is a triangular 
lattice of circular air holes with radius r = 0.3 a and 
lattice constant a^. — a in a dielectric background with 
n = 3. We calculate the band structure for light polarized 
with the H field out of the PC plane {H^ polarization) at 
frequencies a/AG (0, 0.5) in the T — M direction, i.e., at 
normal incidence. Using COMSOL, we calculate the field 
in an 8 period slab of the PC, and we apply our method 



to find the largest three Bloch factors. 



w varies: 



it is 



less than 10~ at low frequencies and less than 10^ at 
high frequencies. 

Fig. [2] summarizes the propagation properties of the 
two/three most dominant modes. The moduli of the 
Bloch factors which quantify how the modes' ampli- 
tudes vary with propagation, are shown in Figs. 2(a) and 
2(b). Below the Wood anomaly, an inspection of A and 
7 shows that the third mode is barely excited by the nor- 
mally incident plane wave, and this reduces the accuracy 
of the results (Fig. 2(a)). Ignoring the uncoupled mode at 
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FIG. 2: (Color online) Complex band structure for the 
PC. The Wood anomaly (a/A = 0.333) is marked. The 
modes are sorted into colors by where two modes 
are propagating (i.e., have |^| = 1), they are sorted by 
|arg(/i)|. (a) Magnitude of Bloch factors with three 
Bloch modes found at all frequencies, (b) with two 
Bloch modes found below the Wood anomaly, three 
above, (c) Argument of Bloch factors, (d) Complex 
band structure in 3D. 



low frequencies (where the p = 1 grating order is evanes- 
cent and so may not be used to excite the structure, 
mentioned in Sec. Ill I 



as 



increases the accuracy of the 
other two modes (Fig. 2(b)). The complex arguments of 
the Bloch factors, which quantify how phase is acquired 
through propagation, are shown in Fig. 2(c), and the in- 
formation about amplitude and phase is summarized in 
a single plot in Fig. 2(d). Aside from slight errors in the 
phase of strongly evanescent modes in Fig. 2(c), there is 
good agreement between Fig. |2] and Bloch factors calcu- 
lated by highly accurate multipole techniques. 

Figure [2] shows that at frequencies below the Wood 
anomaly there is at most one propagating Bloch mode, 
which becomes evanescent in the first bandgap with a 
decay factor \\x\ of no less than 0.5; it still decays far 
more slowly than the other evanescent Bloch modes at 
that frequency. Fig. 2(c) shows that for the evanescent 
modes, either or tt phase is acquired across each unit 
ceU. 



B. Antireflection coating 

Our next application is to reproduce the design of an 
antirefiection coating we presented previously^, found 
using PC impedances calculated with a specialized 
transfer-matrix methocPH. As in this previous paper, our 
design strategy is to try out a very large number of po- 
tential coatings, and choose the coating that gives the 
lowest reflectance off the coated structure. The use of PC 
impedances makes this a feasible problem, as the evalu- 
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ation of each coating is quick, involving a few operations 
on M X M (here 3x3) matrices. 

The target PC is a triangular lattice with lattice con- 
stant ax = a, consisting of air holes in a dielectric back- 
ground with n = 2.86. The holes are cylinders with ra- 
dius r — 0.25 a. We seek to coat the PC to minimize 
reflection for light with frequency a/A = 0.38, incident 
from air at an angle of 30° in the Ez polarization. At 
this frequency and incident angle, Mmin = 2; we consider 
a total of 3 modes to ensure accuracy. As in our previous 
worl?^^, we seek a two-layer coating, where the degree 
of freedom is ay, the lattice vector component perpen- 
dicular to the air/PC interface. For a regular triangular 
lattice, tty = ^a. 

We choose 121 candidate PCs with ay e [0.6, 1.8] ^a 
and simulate 8 periods of each in COMSOL. We apply 
our method to the resulting data, using the Bloch factors 
of the previous PC as the starting estimate for the next. 
BlochCode processes the 121 PCs in approximately 13 
minutes on a 3.06 GHz Intel Core 2 Duo desktop com- 
puter. An equivalent approach that only requires one 
PC to be evaluated is detailed in Sec. |IVC[ we do not 
use it here since the purpose of this section is to demon- 
strate the reliability and consistency of the optimization 
procedure. 

We then calculate the reflectances off the 121^ = 14641 
coated stacks (Fig. |3]) , which takes 34 seconds on a sin- 
gle core of the desktop computer. The optimal coat- 
ing is found to have thicknesses ayi ~ 1.53 ^a and 
a.y2 = 0.65 and reduces the reflectance of the struc- 
ture from R = 0.945 to i? = 1.96 x lO""*. The results in 
Fig. |3] agree well with data calculated by a highly accu- 
rate multipole scattering matrix method: the RMS dif- 
ference is 3.4 X 10"'^, and the only noticeable differences 
occur on the two sharp resonant features near the lower 
edge of the flgure. Specifically, the multipole-based calcu- 
lations show that the coating reduces the PC's reflectance 
from R = 0.943 to i? = 4.29 x 10"*. 



C. All-polarization antireflection coating 

Finally, we apply our methods to find an all- 
polarization antireflection coating for a silicon-based self- 
coUimating square-lattice photonic crystal presented by 
Park et al.^^. They investigated this class of structures 
using a scalar treatment of reflections, and were able to 
design an all-polarization coating at a/A — 0.28, be- 
low the first Wood anomaly. Since their scalar treat- 
ment does not support multiple propagating or evanes- 
cent Bloch modes, it generally does not work above the 
Wood anomaly. Our method does not have this limita- 
tion and we demonstrate this by designing an antireflec- 
tion coating for both polarizations at a frequency well 
above the Wood anomaly, using more than one Bloch 
mode. 

Park et al^ showed that at a/A = 0.368, a 2D sili- 




FIG. 3: (Color online) Reflectance of the coated PC as 
a function of a^i and ay2, the relative thicknesses of the 

two coating layers, calculated using PC impedances 
from BlochCode. The minimum reflectance is marked. 



con {n — 3.518) PC with r — 0.45 a is self-coUimating 
for both polarizations at normal incidence. The large ra- 
dius is an extreme case that is challenging to simulate 
accurately. At this frequency Mmin — 3, so for polar- 
ized light we include AI — 3 modes in our calculations, 
with light incident from the p ~ I grating order so that 
the otherwise uncoupled mode is excited. For light, 
this procedure does not yield accurate results — Bloch fac- 
tors are calculated accurately, but the calculated reflec- 
tion coefficients differ from those calculated directly in 
COMSOL. The calculated impedances prove sufficiently 
accurate to design an effective antireflection coating, but 
the inaccuracies mean that the coating is not optimal. 

To avoid these inaccuracies in Hz polarization, we ex- 
ploit the symmetry that causes the uncoupled mode. The 
physical structure and normally incident field are both 
symmetric about the y-axis, and so modes without even 
symmetry are not coupled to. Therefore we formally ig- 
nore the uncoupled odd mode, in each PC and in the 
reference medium, setting M = 2. In our COMSOL 
simulations for this structure, light is normally incident. 

In Fig. 2 of Park et al.'s papei'^H, they state that 
R ~ 0.28 for Ez polarized light, and R ~ 0.35 for Hz 
light. We calculate with BlochCode that a semi-infinite 
slab of the PC has R = 0.284 for Ez, and R = 0.354 for 
Hz polarized light at this frequency, when incident from 
silicon. Specialized FEM-based transfer-matrix calcula- 
tions agree, showing R — 0.284 for Ez polarization, and 
R = 0.357 for Hz polarization. 

At a/A = 0.368, normally incident light is reflected by 
the PC into three propagating diffraction orders. Due to 
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the symmetries of the problem, the ±1 orders are only 
excited in an even superposition, so light is reflected into 
two modes. A successful coating needs to suppress reflec- 
tion into both these modes simultaneously, and so must 
balance two modes' amplitudes and two modes' phases 
simultaneously for each polarization. Thus the design of 
a perfect all-polarization coating requires 8 continuous 
degrees of freedom. Rather than trying to search an 8- 
dimensional parameter space, which is computationally 
expensive even when the evaluation of each point is effi- 
cient , we consider coatings with four degrees of freedom 
and accept that we are unlikely to find an all-polarization 
coating with zero reflectance. Nevertheless, this is a par- 
ticularly difBcult problem: not only do we need many 
degrees of freedom to find a satisfactory coating, but if 
either of the Bloch factors in a PC is incorrect or any el- 
ement of the PC's impedance matrix is wrong, then the 
calculated net reflection off the structure is incorrect as 
well. 

To limit the coating's thickness, we embed the four de- 
grees of freedom into two rows of holes by varying both 
the hole radii, ri and r2, and the space after the layers, di 
and d2 (Fig. |4]) . Increasing di and ^2 is similar to increas- 



ing ay, as in Sec. IV B but because the candidate PCs 
are independent of d, only one PC per radius needs to be 
simulated in COMSOL. Furthermore, the properties of 
the layers of silicon with thickness di may be calculated 
analytically. We consider 36 possible hole radii in the 
range e [0.10,0.45] a and 99 values of di G (0,1) a. 
To allow a thin coating, we set ay = 2r + 0.1a for each 
PC. If necessary, additional degrees of freedom could be 
added to find a coating with even lower reflectances. 



the large number of candidate coatings 1.3 x 10^), the 
embarrassingly parallel problem was split over 16 cores 
of the workstation, taking approximately 80 minutes per 
polarization. 

The best Ez coating reduces R from 0.284 to 9.56 x 
10~^, and the best Hz coating reduces R from 0.354 to 
3.33 X 10"'*. The best all-round coating is taken to be 
the one with the lowest total reflection in the two po- 
larizations. This coating has ri = 0.13 a, di ~ 0.89 a, 
r2 = 0.17 a, and ^2 = 0.90 a (Fig. In Ez it reduces R 
to 0.0141, and in Hz it reduces R to 0.0197. Calculations 
from a specialized transfer matrix methocP^ agree with 
these results, giving R — 0.0142 in Ez polarization and 
R = 0.0211 in Hz. 

To verify these results without the aid of our special- 
ized methods, implementations of which are not publicly 
available, we simulate the structure using COMSOL Mul- 
tiphysics. Since COMSOL cannot directly calculate re- 
flection coefficients off semi-infinite PCs, we simulate a 
20-period section of the uncoated PC surrounded by the 
background dielectric, and compare the results to a simu- 
lation with the antireflection coating on both sides of the 
PC section. BlochCode calculates the reflectance of the 
uncoated and coated structures to be 0.407 and 0.0124 
respectively in the Ez polarization, and 0.574 and 0.0074 
in the Hz polarization. The COMSOL simulations agree 
with these results, showing that the coating reduces R 
from 0.407 to 0.0129 in the Ez polarization, and from 
0.585 to 0.0055 in the Hz polarization. 



V. DISCUSSION 8i CONCLUSION 
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FIG. 4: Schematic of the all-polarization antireflection 
coating, ri and r2 are the radii of the holes in the first 
two layers, and di and d2 are the thicknesses of the 
extra silicon background layers between the first few 
rows of holes. For this coating, ri = 0.13 a, di ~ 0.89 a, 
r2 ~ 0.17 a, and d2 = 0.9 a. 

On a single core of a 16 x 2.4 GHz Intel Xeon-Quad 
workstation, it took a total of 15 minutes to find the 
modes of the 36 PCs in the two polarizations. For Ez 
polarization, w'^ ~ 10~^ for most radii, and for Hz polar- 
ization ranged roughly from 3x 10"'^ for thin unit cells 
to 10"'' for the thicker cells with larger radius. Due to 



We have detailed a method for calculating the com- 
plex band structure and impedance of PCs. The method 
takes into account structural symmetries in the PC, and 
enforces relationships between the fields of forward and 
backward modes, thus improving the method's accuracy 
by eliminating ill-conditioning and constraining modal 
fields. We have applied the method to three cases, and 
have demonstrated that it works for a variety of square 
and triangular lattice 2D photonic crystals, for light in 
both polarizations and at different incident angles. We 
have demonstrated that our method works at frequencies 
both above and below the first Wood anomaly, the fre- 
quency above which scalar methods cannot adequately 
describe light propagation and reflection in PCs. 

The stronger the excitation of a Bloch mode, the more 
accurately our method calculates its properties. Thus the 
method is well-suited to calculating reflection and trans- 
mission through arbitrary PC stacks, where the most im- 
portant modes are those that are strongly excited. Since 
PC impedances make it so easy to calculate the reflec- 
tion and transmission properties of many combinations 
of PCs in a stack, it is feasible to search large parameter 
spaces of PC stacks for particular reflective properties 
over a range of frequencies, incident angles and polariza- 
tions. The method can be used to design not only all- 



9 



polarization antireflection coatings, but also broadband 
antireflection coating&i-^, polarization filters, angular fil- 
ters, and other devices. 

Ha et al. have applied their method to slab PC 
waveguidefi^. We have not yet applied our method to 
any 3D structure. As long as the x — z plane mirror 
symmetry is present, our method for finding the com- 
plex band structure remains valid. The field of a slab 
waveguide might be sampled only over the PC's surface 
(as in a SNOM experimenlfiSl) or throughout the entire 
volume of the structure (as in a simulation); cither case 
provides sufficient information to determine the modal 
fields within the sampled region and the associated com- 
plex band structure. However, the impedance formalism 
is yet to be developed for 3D structures. 

Our method is also valid for finding modes of PC wave- 
guides, using supercells. Calculation of reflection and 
transmission matrices between PC waveguides is yet to 
be demonstrated using impedances, but they have previ- 
ously been calculated directly from the supercell's E and 
H matrice^. 

Bloch mode analysis is a valuable tool in understand- 
ing light's interactions with PCs. Using an EM solver 
and our method, for which source code is available^, it 
is straightforward to find a PC's complex band structure 
and its impedance. Respectively, these quantities dic- 
tate how the Bloch modes travel through the PC, and 
which modes they couple with at a PC interface. If these 
quantities are known for a set of PCs, then it is fast and 
efficient to calculate how light travels through arbitrary 
stacks of the PCs. 
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